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The stellar mass function of the most massive galaxies at 
3 ^ z < 5 in the UKIDSS Ultra Deep Survey 
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ABSTRACT 

We have analysed a sample of 1292 4.5 fim -selected galaxies at z ^ 3, over 0.6 deg 2 
of the UKIRT Infrared Deep Survey (UKIDSS) Ultra Deep Survey (UDS). Using 
photometry from the U band through 4.5 /im , we have obtained photometric redshifts 
and derived stellar masses for our sources. Only two of our galaxies potentially lie at 
z > 5. We have studied the galaxy stellar mass function at 3 ^5 z < 5, based on the 
1213 galaxies in our catalogue with [4.5] ^ 24.0. We find that: i) the number density 
of M > 10 11 Mq galaxies increased by a factor > 10 between z = 5 and 3, indicating 
that the assembly rate of these galaxies proceeded > 20 times faster at these redshifts 
than at < z < 2; ii) the Schechter function slope a is significantly steeper than that 
displayed by the local stellar mass function, which is both a consequence of the steeper 
faint end and the absence of a pure exponential decline at the high-mass end; iii) the 
evolution of the comoving stellar mass density from z = to 5 can be modelled as 
log 10 p M = -(0.05 ± 0.09) z 2 - (0.22 =F 0.32) z + 8.69. At 3 z < 4, more than 30% of 
the M > 10 11 Mq galaxies would be missed by optical surveys with R < 27 or z < 26. 
Thus, our study demonstrates the importance of deep mid-IR surveys over large areas 
to perform a complete census of massive galaxies at high z and trace the early stages 
of massive galaxy assembly. 
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1 INTRODUCTION 

The study of massive (M > 5 x 10 10 M ) galaxies at high 
(z > 3) redshifts allows for the investigation of the first 
epochs of efficient stellar mass assembly, when the Universe 
was less that a few gigayears (Gyr) old. It is now accepted 
that around 20-40% of the massive galaxies we know to- 
day were already in place by z ~ 2 (Fontana et al. 2004; 
Caputi et al. 2005, 2006a; Daddi et al. 2005; Labbe et 
al. 2005; Saracco et al. 2005; Papovich et al. 2006; Arnouts 
et al. 2007; Pozzetti et al. 2007; Wuyts et al. 2009; Ilbert 
et al. 2010). Above redshift z ~ 3, massive galaxies are 
more difficult to find (e.g. McLure et al. 2006; Kodama et 
al. 2007; Rodighiero et al. 2007) and they become very rare 
by z ~ 5 — 6 (e.g. Dunlop, Cirasuolo & McLure 2007). 

The galaxy populations discovered at z > 6 seem to al- 
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most exclusively consist of intermediate or low mass galaxies 
(M < 10 10 Mq), and their current study is mainly focused 
on constraining the epoch and sources of reionisation (e.g. 
Bunker et al. 2004; Bouwens et al. 2008; McLure et al. 2009, 
2010; Oesch et al. 2010). In fact, with the current instrumen- 
tation, only the rest-frame ultra-violet (UV) emission can 
be observed for the vast majority of these galaxies, which 
is used to estimate their levels of star formation and abil- 
ity to produce and liberate sufficient Lyman-a photons to 
ionise the intergalactic medium. Any constraints on the stel- 
lar masses of z > 6 galaxies are still poor, due to the lack of 
sensitive data at wavelengths that map the galaxy rest-frame 
near-infrared (IR) at these redshifts (although see e.g. Labbe 
et al. 2010 for an attempt). Nevertheless, different pieces of 
observational evidence suggest that massive galaxies as a 
significant population only appear at later times. 

Cosmological models of galaxy formation predict that 
massive galaxies can be quickly formed at high-z in the 
high-density fluctuations of the matter density field (Cole 
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& Kaiser 1989; Mo & White 1996). Thus, determining the 
first epoch of appearance and subsequent rise in the num- 
ber density of massive galaxies with redshift constitutes a 
very important constraint on galaxy formation models. In- 
vestigating the period elapsed between redshifts z ~ 3 and 
6-7 is fundamental for this purpose, as it connects the ear- 
liest stages of galaxy formation after the epoch of reionisa- 
tion, with the better-studied period of galaxy evolution at 
1 < z < 3, where a substantial population of galaxies are 
already massive and host very intense star formation and 
quasar activity. 

A global picture of the evolution of galaxy formation 
and growth can be obtained through the study of the galaxy 
stellar mass function at different redshifts. At low redshifts, 
the shape of the galaxy mass function is known down to 
low mass limits (M ~ 10 s M ; e.g. Cole et al. 2001; Baldry, 
Glazebrook & Driver 2008) and is well-fitted by a double 
Schechter (1976) function (Baldry et al. 2008; Bolzonella et 
al. 2010; Pozzetti et al. 2010). Peng et al. (2010) proposed 
that the Schechter-function shape observed for the galaxy 
stellar mass function up to redshift z ~ 2 can be explained 
as a consequence of mass-driven star-formation quenching 
proceeding proportionally to the galaxy star-formation rate 
(SFR) in M > M* galaxies. This mechanism could be at 
play since earlier times, as preliminary determinations over 
small areas of the sky indicate that the Schechter functional 
form could be suitable to describe the bright end of the 
galaxy stellar mass function up to z ~ 3.5 (e.g. Fontana et 
al. 2006; Kajisawa et al. 2009). 

Selecting galaxies by their rest-frame near-IR light con- 
stitutes a good proxy for a stellar mass selection. Rest-frame 
near-IR wavelengths are relatively unaffected by dust, and 
the corresponding mass-to-light ratios have much smaller 
variations with galaxy age than at shorter wavelengths. For 
galaxies at z 3, the rest-frame near-IR light is shifted into 
observed mid-IR wavelengths. The Infrared Array Camera 
(IRAC; Fazio et al. 2004) on board the Spitzer Space Tele- 
scope (Werner et al. 2004) is currently the most suitable in- 
strument to conduct a mass-selected galaxy survey at high 
redshifts. IRAC has operated at wavelengths 3.6, 4.5, 5.8 
and 8.0 fim until the end of the Spitzer cryogenic mission, 
and still operates in the two shortest-wavelength channels 
in the on-going warm campaign. Spitzer data constitute the 
last opportunity to conduct such mid-IR galaxy surveys un- 
til the advent of the James Webb Space Telescope (JWST) 
after 2014. 

The Spitzer Ultra-Deep Survey (SpUDS; P.I. J. Dun- 
lop) is a Legacy Program that has provided IRAC (Fazio 
et al. 2004) and Multiband Photometer for Spitzer (MIPS; 
Rieke et al. 2004) imaging over more than 1 deg 2 centred on 
the UKIDSS UDS field (P.I. O. Almaini). The UDS is one 
of five on-going surveys which comprise the UKIDSS, and is 
characterised by the existence of deep UV through if-band 
ground-based photometric data over an overlapping area of 
~ 0.60 deg 2 . The availability of deep and homogeneous- 
quality multi-wavelength data in this region makes of the 
UKIDSS UDS one of the most suitable surveys to perform 
galaxy evolution studies at high z, as demonstrated by sev- 
eral works (e.g. McLure et al. 2006, 2009; Dunne et al. 2009; 
Cirasuolo et al. 2007, 2010; Hartley et al. 2010; Ono et 
al. 2010). 

In this work, we present the results of an IRAC 



4.5 fj,m selection and multi- wavelength analysis of galaxies at 
z > 3 in the UKIDSS UDS field. In Section jJU we describe 
the different analysed datasets. In Section A3. II we explain 
in detail our sample selection based on the photometric red- 
shifts of the entire 4.5 /im galaxy catalogue. In Section ^4] we 
compute the galaxy stellar mass function and correspond- 
ing stellar mass densities at 3 ^ z < 5. Later, in Section Sj5] 
we discuss the reliability of our IRAC-selected galaxy can- 
didates at 2 > 5. In Section Sj6] we quantify the fraction of 
our 3 ^ 2 < 5 that would be missed by typical deep optical 
selections. Finally, in Section ^7] we summarise our findings 
and present our conclusions. We adopt throughout a cosmol- 
ogy with H = 70kms~ 1 Mpc" 1 , Q M = 0.3 and tt A = 0.7. 
All quoted magnitudes and colours are total and refer to the 
AB system (Oke & Gunn 1983), unless otherwise stated. 



2 DATASETS 

The UKIDSS (Lawrence et al. 2007) UDS has been con- 
ducted over a field centred on RA=02: 17:48 and DEC=- 
05:05:57 (J2000), and benefits from a range of multi- 
wavelength data from X-rays to radio wavelengths. The field 
is defined by its UKIRT Wide Field Camera (WFCAM) cov- 
erage in the J, H and K bands, which is still in progress. 
The data used in this work corresponds to the fifth data 
release (DR5), and reach 5a depths of 24.0, 23.7 and 23.9 
(2-arcsec-diameter aperture) magnitudes in the J, H and K 
bands, respectively. 

The UDS field has been observed at optical wavelengths 
with SuprimeCam on Subaru (Furusawa et al. 2008). This 
has provided B,V,R,i and 2-band data, with corresponding 
3cr limit (2-arcsec-diameter aperture) magnitudes B — 28.4, 
V = 27.8, R = 27.7, i = 27.7 and z = 26.6. Also, com- 
plementing (7-band data (P.I. O. Almaini) have been ob- 
tained with Megacam on the Canada-France Hawaii Tele- 
scope (CFHT). 

At mid-IR wavelengths, the UDS field has been ob- 
served with IRAC and MIPS for a total time of 124 and 
168 hours, respectively, as part of a Spitzer cycle-4 Legacy 
Program (SpUDS; P.I. J. Dunlop). IRAC observations have 
been carried out with the four IRAC filters, i.e. at 3.6, 4.5, 
5.8 and 8.0 fim. One of the key science drivers of the SpUDS 
program is the possibility of exploiting the combined power 
of optical/near-IR and mid-IR data to study the evolution 
of the galaxy stellar mass function at high redshifts. 

We have used the Spttzer/IRAC 4.5 /im data on the 
UDS field to select a galaxy catalogue at z ^ 3 over a net 
area of 0.60 deg 2 of the UDS field with UV through mid-IR 
coverage. This is the clean area where all datasets overlap, 
and excludes all map edges and regions around bright stars. 
In the next section, we explain in detail the steps of our 
high-2 galaxy selection. 



3 THE SpUDS IRAC 4.5 Aim-SELECTED 
GALAXY SAMPLE AT z> 3 

3.1 Sample selection and multi- wavelength 
photometry 

We extracted a source catalogue from the SpUDS 
4.5 fim image using the software SEXTRACTOR (Bertin & 
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J-[3.6] (AB) 



Figure 1. The (B — J) versus (J— [3.6]) colour-colour diagram for 
our 52,693 4.5 fim -selected sources. Stars form a clearly separated 
sequence on the left-hand side of the diagram. 

Arnouts 1996) with a 'mexhat' kernel. This type of kernel 
is very efficient in crowded fields, as it facilitates source de- 
blending. Considering only the region overlapping the 3.6 
/im map, and excluding edges and regions around bright 
stars, our 4.5 /im catalogue contains 67,937 sources. 

We measured aperture photometry for all our sources 
and obtained aperture corrections using the curve of flux 
growth for isolated stars in the field. Our derived total 
4.5 fj,m magnitudes -referenced as [4.5] hereafter- correspond 
to measured 4-arcsec-diameter aperture magnitudes cor- 
rected by a constant -0.31 mag. This aperture size is usual 
for IRAC photometry (see e.g. Ilbert et al. 2010), as it con- 
stitutes a good balance between directly measuring most of 
the source encircled energy and minimising contamination 
from close neighbours (the IRAC point-spread function full 
width half maximum is ~1.9 arcsec at 4.5 /im). 

We performed simulations to assess the completeness 
and reliability of our 4.5 /im catalogue. To test complete- 
ness, we used the IRAF task 'gallist' to generate a list of 
50,000 artificial objects following a power-law distribution 
between magnitudes 18 and 26. We then created a set of 
100 mock maps based on the real 4.5 /im image, in each of 
which we have randomly inserted 500 of the artificial ob- 
jects (using 'mkobjects' in IRAF). We then ran SExtractor 
on each of these mock maps with the same configuration file 
used for the real image, and checked the fraction of artifi- 
cial sources recovered as a function of magnitude. Through 
this procedure, we determined that our 4.5 /jm catalogue is 
80% and 50% complete to magnitudes [4.5] =22.4 and 24.0, 
respectively. 

We tested the reliability of our catalogue by repeat- 
ing the source extraction procedure on the negative of 
the 4.5 /im image, and considering the fraction of negative 
sources versus magnitude. At [4.5]=22.4 mag, the percent- 
age of spurious sources is below 0.5%. At fainter magni- 
tudes [4. 5]=23. 5-24.0 mag, this percentage rises to around 
10%. However, after imposing that the 4.5 /im sources have 



a counterpart in the independently extracted A'-band cata- 
logue (see below), the fraction of spurious sources becomes 
negligible even at such faint magnitudes. 

We measured 3.6 /im aperture photometry for all the 
4.5 /im -selected sources running Sextractor in dual-image 
mode. The derived total 3.6 /im magnitudes correspond to 
the measured 4-arcsec-diameter aperture magnitudes cor- 
rected by a constant -0.27 mag (as also determined through 
the curve of flux growth of isolated stars). 

To compile the corresponding UV through near-IR pho- 
tometry for our galaxies, we extracted an independent cat- 
alogue based on the UDS K-bsaid image, and ran SExtrac- 
tor on dual-image mode on the U, B, V, R, i, z, J and H- 
band maps, using the position of the if-band sources. In 
these bands, we obtained total magnitudes from aperture- 
corrected 2-arcsec aperture magnitudes in all cases. All mag- 
nitudes have been corrected for galactic extinction. 

We finally cross-correlated the 4.5 /im catalogue (that 
included 3.6 /im photometry) with the Tf-band catalogue 
(that contained [/-band through K-band photometry), with 
a matching radius r = 1.5 arcsec. The final overlapping area 
of all our datasets is 0.60 deg 2 . Our 4.5 /im catalogue with 
if-band counterparts over this area contains 52,693 sources. 

We note that the depth of the near-IR images matches 
very well the depth of the IRAC data in the UDS. Within the 
clean overlapping area of 0.60 deg 2 , the K-band catalogue 
allows us to identify more than 98% of the 4.5 /im sources 
with [4.5] < 22.4 mag. For the deeper [4.5] < 24.0 mag cat- 
alogue, the percentage of identifications is 92%. Our relia- 
bility tests performed on the 4.5 /im catalogues suggest that 
most of the remaining unidentified sources are likely to be 
spurious IRAC sources. 

We excluded galactic stars from our sample via a colour- 
colour diagram. As discussed by McLure et al. (2009), the 
use of the SExtractor stellarity parameter CLASS-STAR 
alone is not a secure way to segregate stars from high- 
z galaxies when using ground-based data, as some of the 
galaxies are compact and could also have large stellarity 
parameters (CLASS_STAR > 0.8 - 0.9). Instead, colour 
segregation is much more reliable. Fig. [1] shows that stars 
form a separate sequence in the (B — J) versus (J — [3.6]) 
colour-colour diagram. Through this colour diagnostic, we 
determined that 2372 out of our 52,693 sources are galactic 
stars. Note, however, that this colour-colour diagram can- 
not segregate red dwarf stars, which are a potential source 
of contamination for high-z galaxy samples (cf. Section ij5}. 

Basically all of the 2372 colour-segregated objects have 
CLASSJ3TAR > 0.8, but they constitute less than a half of 
the total number of sources with CLASSJ3TAR > 0.8 within 
our sample (in our case, we measured the CLASS_STAR 
parameter on the if -band images). After the star separation, 
the final catalogue that we considered for further analysis 
contains 50,321 galaxies. 

3.2 Photometric redshifts and selection of the 
z ^ 3 sample 

We derived photometric redshifts z p hot for our 50,321 
4.5 /im -selected galaxies making use of the photometric data 
in 11 broad bands (U through 4.5pm), and a code largely 
based on HYPERZ (Bolzonella, Miralles & Pello 2000; see 
Cirasuolo et al. 2010 for details). The spectral energy distri- 
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Figure 2. Photometric versus spectroscopic redshifts for 733 Figure 3. Estimated stellar masses versus redshifts for the 1292 

galaxies in our 4.5 fim -selected sample. The overall quality of our galaxies in our z ^ 3 sample. A Salpeter (1955) IMF over stellar 

photometric estimates is <T[|zph t — Zs P ec|/(l+Zs P ec)] = 0.05, after masses M = (0.1 — 100) Mq is assumed, 
excluding catastrophic outliers (3% of the data points). 



bution (SED) fitting has been done using the 2007 version of 
the Bruzual & Chariot templates (Bruzual & Chariot 2003; 
Bruzual 2007), and applying the Calzetti et al. (2000) red- 
dening law with 0.0 ^ Ay ^ 3.0 to account for internal 
extinction. 

As a diagnostic of the quality of our derived z p hot, we 
compared these values with the corresponding real redshifts 
for a subset of 733 galaxies for which spectroscopic data are 
available (Fig. [2). The overall quality of our photometric red- 
shifts is very good: the mean of the (z p hot — z spcc )/(l + z spcc ) 
distribution is -0.01 and the dispersion is a — 0.05 after re- 
moving only 3% of outliers. The final photometric redshifts 
considered for this diagram incorporate all the further tests 
performed on the photometric redshifts that are explained 
below. 

Note that, although unfortunately we do not have spec- 
troscopic confirmation for most of the z ^ 3 galaxies that we 
study in this work, the z p hot versus z spC c diagram gives us 
confidence that our discrimination of lower-redshift galaxies 
is reliable. 

We used the z p hot — z spC c diagnostics to iteratively cor- 
rect the zero-points of our photometric data in different 
bands. The final systematic zero-point corrections we ap- 
plied to our photometry is -0.25 mag in the U band, and 
-0.10 mag for the J, H and K bands. 

Our photometric redshifts are the first criterion to se- 
lect the z 3 galaxy sample. Based solely on the primary 
solutions (i.e. the photometric redshift corresponding to the 
minimum \ 2 value for each galaxy), our sample contains 
1608 z 3 candidates. However, no galaxy should have 
any significant flux blueward of the Lyman-a limit at rest 
A r cst = 912 A. So, as a further test, we checked this condition 
on our z 3 candidates. 

At z > 3.6, the Lyman-a limit is redward of the U- 
band filter. Thus, for accepting a galaxy candidate at z > 
3.6, we imposed that it should be a £/-band dropout (i.e. 
it has a £/-band magnitude below the 2a detection level of 



our catalogue). Following the same criterion, for accepting a 
galaxy at z > 4.5, we imposed that it should be undetected 
in both the U and B bands. After these restrictions, a total 
of 1544 galaxies remain in our z ^ 3 sample. 

On the other hand, the SED-fitting procedure yields 
the x 2 distribution as a function of redshift for each galaxy, 
and with this we can assess the probability P(z) that each 
of the remaining candidates is truly at z ^ 3. For 244 of 
them, the SED fitting indicates that secondary, lower red- 
shift (z scc < 3) solutions cannot be discarded within 3a 
confidence (i.e. Xscc. — X P rim < 9)- Comparison of photomet- 
ric with spectroscopic redshifts for a subset of these objects 
indicates that the secondary, lower-redshift solution is the 
correct one in these cases (in the plot shown in Fig. [21 the 
percentage of catastrophic outliers would more than dou- 
ble if we considered the primary photometric redshifts for 
these objects). Therefore, we decided to exclude from our 
final sample these 244 galaxies for which secondary z < 3 
photometric redshift solutions cannot be discarded within 
3a confidence. This is a conservative criterion that is neces- 
sary to guarantee that we do our science analysis on a secure 
high-z sample. 

The 244 rejected sources constitute ~ 15% of our initial 
z ^ 3 candidate sample. We analyse the potential impact 
that these sources would have in the galaxy stellar mass 
function at 3 ^ z < 5 in Appendix A. 

Our remaining z ^ 3 sample contains 1300 galaxies, 
including ten sources that are z > 5 candidates. We sepa- 
rately discuss these ten sources in Section §5\ but we antici- 
pate that 8 out of these 10 candidates are likely not genuine 
z > 5 galaxies. Our final, secure z ^ 3 sample contains 1292 
galaxies. A total of 1215 out of these 1292 galaxies are within 
the 50% completeness limit of our 4.5 fj,m catalogue, i.e. have 
[4.5] ^ 24.0 mag. 
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4 THE GALAXY STELLAR MASS FUNCTION 

AT 3 sC z < 5 

4.1 The 1/Vmax method and the maximum 
likelihood analysis 

We derived stellar mass estimates for each of our galaxies 
from the best SED fitting performed with the 2007 version 
of the Bruzual & Chariot template library. We used a set of 
templates corresponding to a single stellar population, and a 
collection of exponentially declining star formation histories 
with characteristic times ranging between r = 0.01 and 5 
Gyr, all of them with a solar metallicity. The considered 
age grid was such that the maximum possible age would not 
exceed the age of the Universe at the redshift of each galaxy. 
We convolved each template with the Calzetti et al. (2000) 
reddening law with 0.0 ^ Av ^ 3.0 to account for internal 
extinction. 

The UDS multi-wavelength photometry and, particu- 
larly, the availability of IRAC data that sample the rest- 
frame near-IR light of galaxies up to z > 4, allow us to derive 
reliable stellar mass estimates for all our sources. In all cases, 
the stellar masses are based on the best-fitting SED nor- 
malised between the corresponding K and [4.5] band magni- 
tudes of each galaxy. This produces less uncertainties in the 
stellar mass estimates than allowing the SED to be adjusted 
with the galaxy UV/optical light (cf. Section f4.3p . The de- 
rived stellar masses of our sources versus redshift are shown 
in Fig. [3] We assume a Salpeter (1955) initial mass function 
(IMF) over the stellar mass range M = (0.1 - 100) M Q . 

For the purpose of the stellar mass function calculation, 
we considered only the 1213 galaxies of our sample with 
[4.5] ^ 24.0 mag at 3 ^ z < 5, and divided the sample into 
the following redshift bins: 3.0 ^ z < 3.5, 3.5 ^ z < 4.25 
and 4.25 ^ z < 5.0, which sample comparable comoving 
volumes. 

We first computed the galaxy stellar mass function us- 
ing the 1/V m ax method (Schmidt 1968). Although this tech- 
nique involves data binning, it has the advantage of being 
free of any parameter dependence or model assumptions. Be- 
sides, the normalization is directly obtained from the galaxy 
counts in each stellar mass bin. 

The comoving volume V max considered for each galaxy 
depends on the maximum redshift z max at which the galaxy 
can be observed, given the flux limit of the survey ([4.5] ^ 
24.0 in our case). If this redshift is larger than the max- 
imum redshift of the corresponding redshift bin z m&K ^ 

•Zmax(bin) 5 then Vmax = V m ax(bin) 



V min (bin)- Instead, if 



< z 



max(bin) 3 



then V n 



-V 



min(bin) 



To account 



for the sample incompleteness, we weighed each galaxy by a 
correction factor, depending on its 4.5 /jm magnitude. 

Our results for the stellar mass function calculated with 
the 1 /V m ax method are shown in Fig. [4] (circles) and its val- 
ues given in Table [T] The 80 and 50% strict stellar-mass 
completeness limits, corresponding to a single stellar popu- 
lation formed at z —¥ oo with the catalogue limiting magni- 
tudes [4.5] = 22.4 and 24.0 mag, respectively, are indicated 
with vertical dotted lines in the three panels of Fig. [4] 

The error bars of our stellar mass function are based on 
the result of Monte Carlo simulations and incorporate the 
full probability density distribution P(z) of each galaxy in 
redshift space at z ^ 3. 

We used the P(z) function of each of the 1213 galaxies 



involved in the computation of the stellar mass function to 
create 1000 mock catalogues. Each of these catalogues also 
contains 1213 galaxies with the redshifts randomly assigned 
with a probability given by the P(z) function of each source, 
and constrained by the dropout criteria discussed in Section 
13.21 A new stellar mass is derived so as to be consistent 
with the new random redshift and best-fitting SED in each 
case. We then re-computed the stellar mass function with 
the 1/Vmax method for each mock catalogue, and derived la 
errors for the values originally obtained on the real sample. 
Poisson errors are quoted instead when they are larger than 
the errors derived from the Monte Carlo simulations. 

The comparison of the values given in Table [1] shows a 
significant decrease of the galaxy stellar mass function with 
increasing redshift within the 3 ^ z < 5 interval. For all 
stellar masses, we find that the number of galaxies decreases 
by a factor of ~ 10 between 3.0 ^ z < 3.5 and 4.25 ^ z < 
5.0. This is particularly the case for the three largest stellar 
mass bins, in which our sample is 80% complete at all these 
redshifts. Thus, we can conclude that the decrease in the 
mass-function is a real effect. 

We also performed a second, independent computa- 
tion of the galaxy stellar mass function applying the STY 
(Sandage, Tammann & Yahil 1979) maximum likelihood 
analysis. This is a parametric technique that assumes that 
the shape of the mass function is known. However, in con- 
trast to the 1/Vmax method, this method does not involve 
data binning and does not contain any implicit assumption 
on a uniform spatial distribution of galaxies. The maximum 
likelihood analysis provides a direct calculation of the stel- 
lar mass function, i.e. it does not constitute a fitting proce- 
dure of the stellar mass function obtained with the 1 /V ma x 
method. 

The corresponding maximum likelihood estimator 
reads: 



C[s k \(zi, Mi)i=i,...,if] = 



n 



r 

J la 



<&( Sfe ,M)dlog 10 (M) 



(1) 



Jlog 10 (Af*) 

where the product is made over the i = 1, N galaxies of 
the sample. $(sfc, M) is the adopted functional form for the 
stellar mass function, which depends on the stellar mass M 
and the parameters Sfc. Mq is the minimum stellar mass at 
which the i-th galaxy would be observable, given its redshift 
Zi and the flux limit of the survey. The weighting factors 
Wi allow to take into account sample completeness correc- 
tions. By maximizing C (or, for simplicity, its logarithm), 
one can obtain the values of the parameters yielding the 
maximum likelihood. The normalization factor $* is recov- 
ered after the maximization, by integrating the obtained 
maximum-likelihood mass function without normalization in 
the completeness range of stellar masses of the survey, and 
making it equal to the number density of observed galaxies. 

For the shape of the stellar mass function, we assumed 
the functional form proposed by Schechter (1976): 

3>(M)dlog 10 (M) = $* (il) 1 xexp(-^) dlog 10 (M),(2) 

where a and M* are free parameters. The Schechter func- 
tional form implies that the value of a dominates the shape 
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Figure 4. The galaxy stellar mass function in three redshift bins between z = 3 and 5. Circles show the results of the 1/V ma x method. 
The error bars include Poisson statistics and the effects of the galaxy probability density distributions in redshift space (see text). Vertical 
dotted lines indicate the strict mass completeness -corresponding to a single stellar population formed at z — > oo— that are imposed by 
the 80% and 50% completess limits of our catalogue at [4.5]=22.4 and 24.0 mag, respectively. The thick curve in each panel represents 
the mass function calculated with an STY maximum likelihood analysis, assuming a Schechter function (solid and dashed line, above and 
below the mass completeness limit corresponding to [4.5]=24.0 mag, respectively). The accompanying thin lines indicate the envelope 
of all possible curves within the 1<t confidence limits (see Fig. [5|. The dashed line with the label z ~ corresponds to the local stellar 
mass function computed by Cole et al. (2001). 



Table 1. The galaxy stellar mass function computed with the 1/V ma x method. 
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of the resulting stellar mass function at M < M* . As we 
mentioned in Section iJTJ this functional form well describes 
the shape of the stellar mass function up to at least z ~ 3. 
We briefly discuss the results of exploring other functional 
forms later in this section. 

The thick curves in Fig.[4]show our maximum-likelihood 
stellar mass function obtained with the STY method in each 
redshift bin. The best a and M* parameter values (i.e. those 
yielding the maximum likelihood in each case) are given in 
Table [2l and the 1, 2 and 3<r confidence levels on these pa- 
rameters are plotted in Fig. \5\ These confidence levels as- 
sume that the likelihood function follows a Gaussian distri- 
bution around the maximum, i.e. they are determined by 
Aln(jC) = -0.5, -2.0 and -4.5, respectively. 

The random errors on a, M*, and $* produced by 
the uncertainties in the photometric redshifts -as computed 
with the mock catalogues created following the P(z) dis- 
tribution for each galaxy- are contained within the confi- 
dence levels shown in Fig. [5] The z p hot uncertainties can 
also produce non-negligible systematic effects on the galaxy 
stellar mass function, particularly the so-called Eddington 



bias (Eddington 1913), which we discuss in detail in Section 

EH 

From the results of our maximum likelihood analysis, 
we find the following. On the one hand, the best values of a 
in the three analysed redshift bins are considerably higher 
than the local mass-function slope a = 1.18 ± 0.03 (Cole 
et al. 2001). We observe this steepening of the stellar mass 
function even when the flux limit of our 4.5 /im catalogue 
mainly restricts our analysis to the region around and above 
M* . In fact, this is not surprising: although it is true that the 
value of a determines the slope of the stellar mass function 
in the faint end, it is also important for the bright end, as it 
modulates the exponential decline of the Schechter function 
at M > M*(cf. eq. [2]). Put another way, a flat value a ~ 1 
will only be possible if the stellar mass function decreases 
exponentially at M > M* . 

Our steep a values are closer to those that best de- 
scribe the optical luminosity function at high-z (e.g. McLure 
et al. 2009). Also, Fontana et al. (2006) and Kajisawa et 
al. (2009) found that the slope of the stellar mass function 
increased with redshift up to z ~ 4. 
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Table 2. Maximum likelihood parameter values for the stellar mass function computed with the STY analysis assuming a Schechtcr 
function. The last column lists the values of the comoving stellar mass density obtained by integrating M 5>(M) dlog 10 (M) over stellar 
masses log 10 (Af) = 8.0 to 13. 



Redshift 


a 
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(2.36+°^) x 10 6 



At 3 < z < 4, Fontana et al. (2006) found a « 1.5 
and AT « 8 x 10 10 M©. Note that, although these values 
are lower than those we find in this work, the agreement 
between our stellar mass function and that of Fontana et 
al. (2006) is very good (cf. Section 34.4)1 . This is because the 
values of the parameters a and M* are coupled, and dif- 
ferent combinations of the two can produce a similar shape 
for the stellar mass function around and above M ~ M* 
(where current observations actually constrain the stellar 
mass function). 

At redshifts 2.5 < z < 3.5, Kajisawa et al. (2009) 
obtained a ~ 1.75, a similar value to that we find here. 
Their corresponding characteristic stellar mass is M* w 
2.6 x 10 11 M Q . 

On the other hand, we find that our maximum likeli- 
hood values for M* and $* decrease with redshift within the 
range 3 ^ z < 5, in consistency with our l/V max -method re- 
sults. The comparison of M* and $* is straightforward in 
the 3.0 ^ z < 3.5 and 4.25 ^ z < 5.0 redshift bins, as the 
a values are very similar. Between these two redshift bins, 
M* decreases to less than a half, and <£>* becomes around 
three times smaller. Therefore, our results indicate a fast 
evolution of the stellar mass function in the short (~1 Gyr) 
elapsed time between redshifts z = 3 and 5. 

As an alternative to the Schechter function, we tried a 
double-exponential form in the maximum likelihood analy- 
sis. The two functional forms can only be distinguished in 
practice when the stellar mass function is well constrained 
at the very bright end (M 3> M*\ cf. Caputi et al. 2007 for a 
discussion in the case of luminosity functions) . In our present 
analysis, we find that the maximum likelihoods obtained 
with either the Schechter or the double-exponential func- 
tions are consistent within 3<r (the resulting stellar-mass- 



function curves are almost identical within the stellar mass 
range sampled in our catalogue). 

A single power-law can be discarded as the shape of the 
stellar mass function at 3.0 ^ z < 3.5 and 3.5 ^ z < 4.25 
with > 6a and > 3<j confidence, respectively. Instead, we 
cannot exclude the single power-law shape at 4.25 ^ z < 5.0, 
for which the maximum likelihood is within ~ 2a confidence 
of the Schechter function case. It is likely that this effect is 
produced simply because we might not be well sampling the 
M < M* regime of the stellar mass function at these red- 
shifts. However, this could alternatively be suggesting that 
the shape of the stellar mass function is actually changing 
around this cosmic epoch. Deeper mid-IR observations that 
allow a better sampling of the faint end of the stellar mass 
function are necessary to confirm that this is a real effect. 

4.2 Tests on the faint-end slope a 

We performed a few tests on our galaxy stellar mass function 
to assess the robustness of the large values we get for the 
faint-end slope parameter a. 

Firstly, we simply repeated our ML analysis fixing the 
value of a to different smaller values. Fig. [5] compares the 
galaxy stellar mass function resulting from alternatively con- 
sidering fixed ot — 1.3, 1.6, and leaving a as a free parameter 
(cf. Section ij4.ip . As before, the normalisation $* in each 
case has been obtained a posteriori, by imposing that the 
integral of the stellar mass function is equal to the source 
number density down to the mass completeness limit. 

Our results show that, at 3.0 ^ z < 3.5 and 3.5 ^ 
2 < 4.25, the best galaxy stellar mass function that we 
obtained with a Schechter function with fixed a = 1.6 is 
still in a reasonable agreement with the stellar mass func- 
tion computed with the 1/V max method. Instead, a fixed 
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Figure 6. Comparison of the galaxy stellar mass function ob- 
tained with the ML analysis at different redshifts, assuming a 
Schechter function with fixed slope a = 1.3, 1.6 and a free a 
value. Circles show the results of the 1/V max method (cf. Fig.[4]|. 



a = 1.3 value significantly under-predicts the stellar mass 
function high-mass end, and at the same time over-predicts 
the number density of galaxies at intermediate stellar masses 
M ~ 2 - 6 x 1O 1O M . These results indicate that a slope 
as small as a = 1.3 is not adequate to describe the galaxy 
stellar mass function derived from our datasets. 

At 4.25 ^ z < 5.0, the same effect is present, although 
the constraints on the stellar mass function provided by our 
data in this redshift bin are less tight, and the case of a small 



Figure 7. The galaxy stellar mass function at 3.0 ^ z < 3.5 
obtained by applying the maximum likelihood analysis only to 
those galaxies with [4.5] ^ 22.4 mag. Circles show the results of 
the 1/Vmax method for all galaxies down to [4.5] = 24.0 mag. 



a value cannot be rejected. (See that the value a = 1.3 is 
contained within the 2a confidence level shown in the right- 
hand panel of Fig[5} • 

Another test we performed to assess the robustness of a 
large a value consisted in restricting our ML analysis to our 
galaxy sample with 80% completeness, i.e. only to sources 
with [4.5] ^ 22.4 mag. We only performed this test in the 
3.0 ^ z < 3.5 redshift bin, where we have a sufficient number 
of [4.5] ^ 22.4 galaxies as to allow both a and M* to be 
treated as free parameters. 

Figure [7] shows our results. Our maximum likelihood 
analysis applied only to the [4.5] ^ 22.4 mag galaxies yields 
a significantly larger value of the faint-end slope, i.e. a f» 2.4, 
and M* « 1.8 x 10 11 M©. However, we see that the extrapo- 
lation of this stellar mass function curve would also overpre- 
dict the 1/Vmax points obtained for stellar masses lower than 
the strict completeness limit imposed by the [4.5] = 22.4 
magnitude cut. 

In conclusion, our tests confirm that the galaxy stel- 
lar mass function at 3 ^ z < 5 is characterised by a high 
value of the slope a. The best slope values we obtain from 
our datasets are 1.8 < a < 2.1, but in practice somewhat 
more modest values such as a — 1.6 could still be adequate. 
Instead, much smaller a values as those characterising the 
galaxy stellar mass function at low-z can be clearly dis- 
carded. 

We note that, in the Schechter function, the a value de- 
termines the slope of the stellar mass function at the faint 
end, but has also the role of modulating the exponential de- 
cline at M > M* . Thus, a large a value in the Schechter 
function is indicating the absence of a pure exponential de- 
cline at the high-mass end. This effect is particularly at play 
within our sample, as it mostly constrains the stellar mass 
function around and above M* rather than M M* . 



4.3 Analysis of the Eddington bias 

In this section, we investigate in more detail the effect that 
the uncertainties in the z p hot and stellar mass estimates 
could produce on the galaxy stellar mass function. These 
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Figure 8. The relative stellar mass error versus the estimated 
stellar mass of each of our galaxies. These errors consider the 
parameter degeneracies of the SED fitting at the best-fitting red- 
shift. In all cases, to derive the stellar masses, the SED model 
is normalised between the K and [4.5] band magnitudes of each 
source. 



uncertainties are the consequence of degeneracies produced 
in the SED fitting, and depend on different sources of errors, 
including the photometric errors, the wavelength coverage, 
and the limited SED template grid (see e.g. Kitzbichler & 
White 2007; Fontanot et al. 2009; Behroozi et al. 2010, for 
further discussions on this subject). 

In particular, we focus here on the analysis of the Ed- 
dington bias (Eddington 1913), which can affect the bright 
end of the galaxy luminosity or stellar mass function, due to 
its exponential decline. In simple terms, the z p hot and ad- 
ditional stellar mass uncertainties produce a scatter in the 
derived stellar mass value of each galaxy. But, even if this 
scatter were symmetric around the best stellar-mass value, 
the effect on the steep bright end of the galaxy stellar mass 
function might not be, as the increment in the number den- 
sity of sources is relatively more important at higher than 
lower stellar masses. As a consequence, this effect can pro- 
duce an artificial flattening of the stellar mass function at 
the high-mass end. 

For simplicity, and to concentrate exclusively on the 
Eddington bias effect, we performed our analysis considering 
that the possible z p hot for each galaxy follow a Gaussian 
distribution around the best z p hot- In each case, the standard 
deviation a is given by the 68% confidence values obtained in 
the SED fitting procedure. We recall that the overall quality 
of our photometric redshifts is cr[|z phot - z spoc |/(l + z spcc )] = 
0.05, after excluding catastrophic outliers (which constitute 
3% of the sample with spectroscopic redshifts). These errors 
in the z p hot directly produce variations in the derived stellar 
masses. However, in addition, we considered a further error 
on the stellar mass values, produced by degeneracies in the 
template parameter space (SED type, age, reddening), at a 
fixed redshift. Fig[8] shows the relative stellar mass error for 
each of our galaxies that is produced by the degeneracies of 
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Figure 9. The effect of considering the z p h t and stellar mass 
errors of each galaxy on the determination of the galaxy stellar 
mass function. In each panel, the thick line corresponds to the 
stellar mass function obtained applying the maximum likelihood 
analysis on the real sample of 1213 z ^ 3 galaxies with the original 
redshifts (cf. Fig- HJ). The multiple thin lines correspond to 100 
Monte Carlo realizations, with the maximum likelihood analysis 
performed on 100 mock catalogues that contain 1213 galaxies 
each. The redshift of each source in the mock catalogues has been 
randomly assigned following a Gaussian distribution around the 
best Zp^ot- The stellar mass corresponds to a random value taken 
from another Gaussian distribution, which is centred at the stellar 
mass given by the random z p hot- 
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the SED template fitting. These errors are < 50% in almost 
all cases, and < 30% for ~90% of the most massive galaxies. 

We remind that we computed the stellar mass of each 
galaxy from the best-fitting SED, which is normalised freely 
between the corresponding K and [4.5] band magnitudes. In 
this way, the variations in the stellar mass estimates due to 
degeneracies in parameter space are minimised. Stellar mass 
estimates which are based on the normalisation of the SED 
adjusted on UV/optical bands carry much larger uncertain- 
ties. In any case, the errors on the stellar mass estimates 
produced by the uncertainties in the z p hot are dominant, 
and these are governed by the quality of the photometry in 
all the different bands. 

We created a set of mock galaxy catalogues with 1213 
galaxies each, where the z p hot and stellar masses have been 
randomly assigned considering two independent Gaussian 
distributions, as explained above. The Gaussian distribution 
on the stellar masses has been centred at the stellar mass 
corresponding to the random z p hot in each case, and has a 
dispersion produced by the template parameter degeneracy 
at a fixed redshift. 

Fig. shows the galaxy stellar mass function obtained 
with the maximum likelihood STY analysis over 100 Monte 
Carlo realizations, in different redshift bins. The original 
stellar mass function obtained with the real z > 3 sample is 
shown for comparison in all cases. 

Both at 3.0 ^ z < 3.5 and 3.5 ^ z < 4.25, we see that 
the galaxy stellar mass function computed on the mock cat- 
alogues quite symetrically scatters around the real one at all 
stellar masses. This shows that the z p hot /stellar mass errors 
do not introduce any significant bias on the determination 
of the galaxy stellar mass function. This is not surprising, 
as the combination of large a and M* values indicates that 
the stellar mass function is characterised by a not so steep 
decline in the high-mass end. We also note that the scatter 
of the galaxy stellar mass function produced by the these 
errors is contained within the la confidence envelope curves 
shown in Fig(4] 

In the redshift bin 4.25 ^ z < 5.0, instead, we see that 
the simulated stellar mass function curves can be flatter than 
the real one at the high-mass end. As the high-mass end is 
constrained by less galaxies at such redshifts, the variations 
in the stellar mass function produced by the z p hot /stellar 
mass errors are much more drastic. This effect is a mani- 
festation of the Eddington bias and, in our case, it appears 
to affect our 4.25 ^ z < 5.0 galaxy stellar mass function at 
M > 1.5 x 10 11 M . 

Correcting the original stellar mass function for the 
Eddington bias requires some assumption on how the true 
function can be inferred from the observed one. Following 
Eddington (1940; see also Teerikorpi 2004), the true stellar 
mass function can be approximated by the observed one, 
convolved by a Gaussian kernel: 



$ rea! (M) 



$ oiM -(M) x e" 



yCT P 



(3) 



where a is the characteristic dispersion observed in the 
stellar mass. The parameter j3 corresponds to the slope of 
the stellar mass function bright end, i.e. we can consider 
/3 ?s -M/M*. 

Taking a = 0.3 (which is roughly the dispersion of 
AM/M* in our case), we get a correction of 0.03 dex 
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Figure 10. Comparison of our stellar mass functions at 3.0 ^ 
z < 3.5 and 3.5 ^ z < 4.25 with those derived by other au- 
thors from near or mid-IR-selected galaxy samples [Fontana et 
al. (2006; F+06); Perez-Gonzalez et al. (2008; PG+08); Kajisawa 
et al. (2009; K+09); Marchesini et al. (2009; M+09) with their 
model set 8, which corresponds to the same template library, 
metallicity and reddening law as those we adopt here]. In the 
upper and lower panels we show the results obtained with the 
1/V m ax method and the maximum likelihood analysis, respec- 
tively. 



for our 4.25 z < 5.0 galaxy stellar mass function at 
M = 1.5 x 10 11 M B , and 0.13 dex at M = 3 X 10 11 M H . 



4.4 Comparison with other works 

It is instructive to compare our stellar mass function with 
those obtained in previous studies from both near and mid- 
IR galaxy-selected samples. This is possible in the redshift 
range 3 < z < 4, where our analysis and other existing 
studies overlap. 

Fontana et al. (2006) and Kajisawa et al. (2009) studied 
the galaxy stellar mass function up to z ~ 3.5—4.0 within the 
Great Observatories Origins Deep Survey South and North 
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fields (GOODS-S and GOODS-N), respectively, using K- 
band galaxy selections. Also, Perez-Gonzalez et al. (2008) 
and Marchesini et al. (2009) performed similar studies on 
larger areas (~ 664 and 511 arcmin 2 , respectively). The re- 
sults of these works are compared with our own in Fig. [TO] 
All stellar masses in this figure refer to a Salpeter (1955) 
IMF over stellar masses M — (0.1 - 100) M . Stellar masses 
referring to this IMF are a factor of 1.6 higher than those 
in the pseudo-Kroupa (2001) IMF adopted by other authors 
(such as e.g. Marchesini et al. 2009). 

From the two panels of this figure we can see that our 
estimation of the galaxy stellar mass function at 3 ^ z < 4 
(calculated over an area >2000 arcmin 2 ) is in very good 
agreement with those obtained by Fontana et al. (2006) and 
Kajisawa et al. (2009), and in marginal agreement with that 
of Marchesini et al. (2009). Note that, for a comparison with 
Marchesini et al. (2009), we are considering their stellar mass 
function derived using the Bruzual & Chariot (2007) tem- 
plates, i.e. the same SED library we adopt in this work. 
We refer the reader to the Marchesini et al. (2009) paper 
for an analysis of the stellar mass function variations pro- 
duced when using different template libraries, metallicities 
and reddening laws. 

Our stellar mass function is also consistent within the 
error bars with that of Perez-Gonzalez et al. (2008) around 
a stellar mass M ~ 10 11 Mq. However, the stellar mass func- 
tion determined by these authors is significantly discrepant 
with that of Fontana et al. (2006) and our own at the high- 
mass end (log 10 M > 11.4). Although these differences might 
be simply due to sample variance effects, we note that our 
surveyed area is around three times larger than that consid- 
ered by Perez-Gonzalez et al. (2008). Thus, we expect that 
our determination of the high-mass end of the stellar mass 
function is more robust. Moreover, our photometric reshifts 
have better precision and smaller fraction of catastrophic 
outliers. This could also explain the discrepancies, because 
the potential inclusion of a few low-redshift contaminants in 
their 3 < z < 4 sample could be resposible for the enhance- 
ment they found at the high-mass end of the stellar mass 
function. 

The comparison of the best parameter values for 
the maximum likelihood analysis performed with different 
datasets should be done with care, as the different free pa- 
rameters are coupled, and different sets of parameter values 
can produce a similar Schechter function. Keeping this in 
mind, it is particularly interesting to note the differences 
found by different authors in the slope a that governs the 
shape of the stellar mass function at the faint end. Perez- 
Gonzalez et al. (2008) determined that the best value of this 
slope at3<z<4isQ~ 1.2, i.e. close to the local value, 
while Fontana et al. (2006) and Kajisawa et al. (2009) ob- 
tained a ~ 1.5 and 1.75, respectively. In this work, we find 
an even higher value for the slope a = 1.86. As we discussed 
in Section £|4.2I smaller values such as a ~ 1.6 would still 
be suitable to describe our galaxy stellar mass function, but 
q < 1.3 are clearly rejected from our data. Unfortunately, a 
more precise constraint on the a value requires a better sam- 
pling of the galaxy stellar mass function faint end and, thus, 
will not be possible until complete, deeper mid-IR surveys 
become available. 



4.5 The evolution of the stellar mass density with 
redshift 

We integrated our resulting stellar mass distribution func- 
tion M $(M) dlog 10 (M), over stellar masses log 10 (M) = 8.0 
to 13.0, in order to obtain estimates of the total comoving 
stellar mass densities (pm) at different redshifts. Our results, 
along with other recent determinations from the literature, 
are shown in Fig. [TTJ 

Our derived comoving stellar mass densities are 
logio(PAr) = 7.48l°;22, 7.05l£;iJ and 6.37±°-l£, at redshifts 
z ~ 3.25, 3.9 and 4.6, respectively (where pu is in units 
of Mq Mpc -3 ). These values correspond to ~6, 2 and 0.5% 
of the local value. Note that, even at the highest redshifts, 
the Eddington bias effect discussed in Section i|4.3| produces 
a negligible correction, as the contribution of galaxies with 
stellar masses M > 1.5 x 10 11 Mq to the stellar mass density 
at z ~ 4.6 is only ~ 4%. 

Our results are in agreement with those obtained by 
Fontana et al. (2006) at 3 < z < 4, Kajisawa et al. (2009) 
at 2.5 < z < 3.5, and Perez-Gonzalez et al. (2008) at 3.0 < 
z < 3.5, within the error bars (note that, for the considered 
a values, the main contribution to the stellar mass density 
integral comes from the region M < M* , and this is why 
Perez-Gonzalez et al. value is consistent with the others in 
spite of the significant differences in the high-mass end of 
their stellar mass function). 

Instead, our comoving stellar mass density at 4.25 ^ 
z < 5.0 is significantly lower than the best estimates of Stark 
et al. (2007) and McLure et al. (2009) at z ~ 5, based on 
the study of deep z-band-selected galaxy samples (although 
we note that our value is still consistent with the lower limit 
of p M = 10 6 M Q Mpc" 3 inferred by Stark et al.). We note 
that, although our galaxy stellar mass function at these red- 
shifts has a similar slope a to that obtained by McLure et 
al. (2009), the values of M* and $* are significantly dif- 
ferent, explaining the differences observed in the integrated 
stellar mass densities. 

The high uncertainty still existing in the value of the 
comoving stellar mass density at z > 4 is due to different 
factors. Current mid-IR surveys such as the one analysed 
here are not sufficiently deep to well sample the stellar mass 
function at M < M* . Instead, optically-selected galaxy sam- 
ples can reach relatively fainter fluxes and, thus, are more 
complete for tracing less massive galaxies at high redshifts. 
However, for most of these less massive galaxies, the rest- 
frame near-IR emission is unknown, with the result that 
the stellar mass estimates obtained in optical surveys are 
very uncertain. Therefore, as we concluded in Section £|4. 1 1 
deeper mid-IR observations will be fundamental to better 
constrain the stellar mass content locked in galaxies when 
the Universe was less than 1.0-1.5 Gyr old. 

Considering our results along with those compiled from 
the literature for different redshifts (see Fig. Hip , we can 
model the evolution of the stellar mass density from z ~ 5 
to the present day. We perform a least-squares fitting for 
the data points, assuming a simple functional form for the 
redshift evolution: 

logio Pm = az2 + bz + c, (4) 

where a and b are free parameters, and c is a constant fixed 
to the local stellar-mass-density value: c = 8.69 (Cole et 
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Figure 11. Total comoving stellar mass densities versus red- 
shifts. The two curves represent the least-squares fittings to the 
data points [dashed and dotted lines, respectively for the cases 
of excluding and including the stellar mass densities estimated at 
z ~ 5 (crosses), which are based on optical surveys]. The ref- 
erences in the plot label are the following: Cole et al. (2001; 
C+01); Fontana et al. (2004, 2006; F+04, F+06); Caputi et 
al. (2006a; C+06); Perez-Gonzalez et al. (2008; PG+08); Kaji- 
sawa et al. (2009; K+09); Stark et al. (2007; S+07); McLure 
et al. (2009; McL+09). All stellar mass densities in this dia- 
gram correspond to a Salpeter (1955) IMF over stellar masses 
M = (0.1 - 100) M . 



al. 2001). The results of our least-squares fitting yields: a — 
-0.02 + 0.07 and b = -0.30 + 0.28, when we consider all the 
data points; and a = -0.05+0.09 and b = -0.22+0.32, when 
we consider only the results of near and mid-IR surveys (i.e. 
all the stellar mass density estimates except those at z ~ 5 
derived by Stark et al. 2007 and McLure et al. 2009). 



5 IRAC 4.5 ^m-SELECTED GALAXY 
CANDIDATES AT z > 5 

As we explained in Section H3.ll our final sample contains 10 
sources whose best z p hot indicate that they should be z > 5 
galaxies. These objects survived all the additional tests we 
have performed to ensure that we only considered a reliable 
high-z sample for our analysis. However, it is known that 
the SEDs of z > 5 galaxies can be well mimicked by that of 
other kinds of objects. The main containants are M, L and 
T dwarves and intermediate-redshift extremely red galaxies 
(ERGs; see e.g. Caputi et al. 2004; Simpson et al. 2006). 
Note that M, L and T dwarves cannot be easily segregated 
by making use of a colour-colour diagram such as that shown 
in Fig. [I] 

As discussed by McLure et al. (2006), an effective dis- 
criminant for the z > 5 galaxy population against lower- 
z galaxies is the colour-colour criterion (R — z) ^ 3 and 
(z — J) ^ 1, but this is not totally effective to separate out 
M dwarf stars. 

Among our ten z > 5 candidates, one is an i-band 
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Figure 12. Postage stamps for the two candidates in our 
4.5 fim sample that are likely genuine z > 5 galaxies: 
uds_ir45.51741 (left) and uds_ir45.96606 (right). The top and 
middle panels show the z-band and 4.5 /xm-band images of these 
sources, respectively. Each stamp corresponds to an area of 10 X 10 
arscec 2 . The bottom panels show the correspoding best-fitting 
SED. 



dropout and another one is a J-band dropout within the 
UDS data depth. These two sources have the typical colours 
of ERGs, i.e. (i-K) > 4 or {R-K) > 5.3 (Vega), and both 
are detected in the Spitzer/MIPS 24 )im images with a flux 
density Sv(24/im) w 160 pJy. Given these properties, we 
conclude that these two galaxies are very likely dusty star- 
bursts at intermediate redshifts rather than z > 5 galaxies. 

Other seven out of the ten z > 5 candidates in our 
sample have colours (R — z) < 3. Six of these are either 
ERGs defined by (i - K) > 4 or (R - K) > 5.3 (Vega), or 
have an SED that is well fitted with the template of an M 
dwarf star (from the Rayner, Cushing & Vacca 2009 library). 

Instead, there is one of these seven sources (id 
udsjr45_51741, with z p hot = 5.07 in our catalogue) that 
is neither an ERG nor any kind of dwarf star. The left-hand 
panels of Fig. [12] show the z-band and 4.5 /an images of this 
source, as well as its best-fitting SED. Although other ob- 
jects present in the field are close to this galaxy, they are 
at a distance d >2 and d > 4 arcsec, in the z-band and 
4.5 /im images, respectively. So, the aperture photometry for 
udsjr45_51741 should not significantly be contaminated by 
the light of neighbour sources. Therefore, we believe that 
udsjr45_51741 could be a genuine z > 5 galaxy in spite of 
its (R — z) < 3 colour. If it were the case, this would be a 
rare example of an old and massive galaxy present in the 
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Figure 13. The percentage of galaxies with [4.5]^ 22.4 (i.e. 
within our 80% completeness limit) at 3 ^ z < 4 that have 
R > 26.6 and z > 25.6 AB total magnitudes. These galaxies 
would be missing in current typical optical surveys, which are 
brighter than these limits. 



early Universe: the best-fitting SED suggests an age of ~ 1 
Gyr and a stellar mass M « 3.8 x 10 11 M Q . 

Finally, one of our 10 z > 5 candidates (id 
udsJr45.96606) has (R - z) = 3.03 ± 0.56 and (z - J) = 
0.51 ± 0.37. The estimated redshift of this source is z p hot = 
5.23. The fitting of its SED with an M dwarf star template 
can be rejected with more than 3<r confidence (with respect 
to the fitting with a galaxy template). This is a well isolated 
object in the optical bands, and does not show any sign of 
blending in the 4.5 /im band (see Fig. 1120. Thus, it is very 
likely another genuine z > 5 source. Our SED fitting indi- 
cates that this is a ~0.2 Gyr old galaxy with a stellar mass 
M « 4.6 x 10 10 M Q . 



6 MID-IR VERSUS OPTICAL SELECTION OF 
HIGH- REDSHIFT GALAXIES 

The bias introduced by doing a flux-limited 4.5 /im selection 
of z 3 galaxies is quite clear: the most massive galaxies will 
be included in the sample, while there is a poor constraint 
on the less massive galaxy populations. Optical selections 
that map the rest-frame UV light of z ^ 3 galaxies, instead, 
favour the selection of sources with high levels of star for- 
mation and/or little dust obscuration, but their biases with 
respect to a stellar mass selection are not obvious. Our aim 
in this section is to use our z ^ 3 sample to investigate this 
problem. 

Current optical surveys can reach faint limits, usually 
R — 27 or z = 26 AB magnitudes (measured through 2- 
arcsec-diameter apertures) over reasonable-size fields, such 
as the UDS or the COSMOS field. So, within our sample, 
we searched for galaxies that would be missed even in these 
deep optical surveys. For a clear comparison, we limited the 
analysis to our 4.5 /jm catalogue with 80% completeness, i.e. 
[4.5] 22.4, and all redshifts 3 ^ z < 4, to be able to explore 
a more or less wide stellar mass range. 



Fig. [13] shows the percentage of our [4.5]^ 22.4 
galaxies at 3 z < 4 that have R > 26.6 and 
2 > 25.6 AB total magnitudes (roughly equivalent to 
R > 27 and z > 26 2-arcsec-diameter aperture mag- 
nitudes), i.e. that will be missed even by deep optical 
selections. The total number of 4.5 /im galaxies analysed 
in each stellar mass bin are 69, 53, 41, 12 and 5 for 
log 10 Af SE [10.7; 10.9), [10.9; 11.1), [11.1; 11.3), [11.3; 11.5) 
and [11.5; 11.7), respectively. 

Our results show that typical deep optical surveys miss 
a significant fraction of massive galaxies at z — 3 — 4. The 
percentage of missed galaxies clearly increases with stellar 
mass: it is (20.3±4.8)% for galaxies with M ~ 6 x 10 10 M , 
and as high as (60 ± 22)% for the rare M ~ 3 - 4 X 10 11 M 
galaxies. 

Note that, although the exact figures depend on the dif- 
ferent magnitude cuts chosen for the R and z bands, our 
conclusions are still valid even when considering slightly 
deeper magnitudes. For example, the fraction of galaxies 
with [4.5] < 22.4 at 3 sC z < 4 that have R > 27 and 
z > 26 AB total magnitudes range from (11.6 ± 3.9)% to 
(29.3±7.0)% and (60±22)%, for stellar masses M ~ 6x 10 10 , 
1.6 x 10 11 and M ~3-4x 10 11 M©, respectively. 

The fact that the fraction of sources missed by deep 
optical surveys increases with stellar mass is directly related 
to an increase in internal extinction: the median extinction 
of our [4.5] < 22.4 galaxies at 3 < z < 4 rises from A v = 0.80 
for stellar masses log 10 M £ [10.7; 10.9), to Av = 2.20 for 
the 5 galaxies with log 10 M £ [11.5; 11.7). 

Our results suggest that, at higher redshifts, when 
galaxies are in general less massive and reddened, deep opti- 
cal surveys should be less biased in selecting the most mas- 
sive galaxy populations. At z < 4, the epoch of obscured 
star formation and black-hole activity is already on-going, 
and this is when IR surveys make a key contribution in dis- 
covering galaxy populations that are usually not detected 
otherwise. 

A total of 4 out of the 5 galaxies with log 10 M € 
[11.5; 11.7) in our 3 ^ z < 4 sample are detected in the 
Spitzer/MIPS 24 /an band with flux densities S v (24 /im) > 
95 /iJy, indicating that these galaxies indeed host dust- 
obscured activity. For two of these galaxies, in particular, the 
24 /im flux densities are quite high (S'„(24/im) > 300 /iJy), 
which at these redshifts could suggest the presence of an 
obscured active galactic nucleus (AGN; see e.g. Desai et 
al. 2008). However, we note that, for none of these two galax- 
ies, the UV through near-IR SED shows any obvious sign of 
an AGN component, i.e. we get good SED fittings using only 
stellar templates (Fig. I14[). 



7 SUMMARY AND CONCLUSIONS 

We have performed a survey of the most massive galax- 
ies present at z ^ 3, over an area of 0.6 deg 2 of the 
UKIDSS UDS field. To have the best possible proxy for a 
stellar mass complete sample, we made our selection in the 
Spitzer/IRAC 4.5 /jm band, which maps rest-frame near-IR 
wavelengths at these high redshifts. We followed up our mas- 
ter 4.5 /^m catalogue of 50,321 sources in 10 broad bands, 
from the [/-band through the IRAC 3.6 /im channel. The 
multi-wavelength follow up has allowed us to model the 
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Figure 14. The best-fitting SED for two galaxies with M > 3.2 X 
10 11 Mg and 5„(24/xm) > 300 /xJy in our 3 < z < 4 sample. The 
galaxy corresponding to the left-hand panel would be detected in 
typical optical surveys, while that of the right-hand panel would 
be missed. 



SEDs of all our galaxies and, with this, obtain redshift es- 
timates and derive stellar masses. Our final sample consists 
of 1292 galaxies at redshifts 3.0 ^ z < 5.23. 

The main goal of our work was the study of the galaxy 
stellar mass function at 3.0 ^ z < 5.0, particularly the evo- 
lution of its high-mass end within this redshift range. Our 
deep and homogeneous datasets over a large field are partic- 
ularly suitable for this purpose. We have found the following: 

• the comoving number density of the most massive 
galaxies (M > 5 x 10 10 and M > 1 x 10 11 M ) declines 
by a factor > 10 within the ~1 Gyr of elapsed time be- 
tween redshifts z — 3 and 5. In turn, the number density of 
M > 10 1 M© galaxies at z — 3 is approximately a factor 
of ten lower than the density ~1 Gyr later at z — 2, and a 
factor of ~ 50 lower than in the local Universe (cf. Caputi 
et al. 2006a). These results altogether indicate that massive 
galaxies have assembled at a very fast rate when the age of 
the Universe was between ~ 1 and 3 Gyr, and this assembly 
has significantly slowed down afterwards. 

Our conclusions on the rapid evolution of massive galaxy 
growth at 2 < z < 5 are fully consistent with the intense star 
formation and black-hole activity characterising this cosmic 
epoch. At redshifts z ~ 2 — 4, a substantial fraction of mas- 
sive galaxies are ultra-luminous infrared galaxies (ULIRGs; 
see e.g. Caputi et al. 2006b; Daddi et al. 2007), and/or host 
AGN (e.g. Alexander et al. 2005; Fiore et al. 2008; Yamada 
et al. 2009). Further investigation of this period appears then 
to be critical to understand the buildup of the most massive 
systems and their evolution to the present day. 

• a Schechter function is still the most suitable form to de- 
scribe the galaxy stellar mass function at z ~ 4. However, at 
4.25 ^ z < 5.0, a single power law cannot be discarded with 
the maximum likelihood analysis. Although this shape de- 
generacy could be an effect produced by the more restricted 
sampling that we have around M ~ M* at such redshifts, 
it could also be an indication that the shape of the stellar 
mass function is changing around that epoch. In fact, the 
transformation of the stellar mass function shape is mani- 
fested by the progressive steepening of the a value that we 
discuss below, and the evolution into a single power law will 
be naturally expected. The possibility of exploring a wider 
range of stellar masses with deeper mid-IR surveys will be 
key to probing the transformation of the galaxy stellar mass 
function shape at high redshifts. 



• the absolute value of the stellar-mass- function faint-end 
slope a increases with redshift. A similar conclusion has been 
previously obtained by Fontana et al. (2006), Kajisawa et 
al. (2009) and Marchesini et al. (2009). In particular, in this 
work we find that the best values are 1.8 < a < 2.1 at 
3 ^ 2 < 5, and that values as small as a < 1.3 can clearly 
be rejected. 

We remark that a steep a value in the Schechter function 
is a consequence of both the stellar mass function low-mass 
and high-mass ends. A flat a ~ 1, as that characterising the 
local stellar mass function, is only possible if the stellar mass 
function has a simple exponential decline at the bright end. 
The rise in the a value observed with redshift is not only 
a consequence of a steepening of the faint end (which, in 
fact, is poorly constrained by all current surveys), but also 
indicates a transformation in the high-mass end of the galaxy 
stellar mass function, with the pure exponential decline no 
longer being observed at high z. 

The steep a values we find at 3 ^ z < 5 are also in 
line with those obtained for the rest-frame UV luminosity 
function at z > 5. We note, however, that the limits of our 
mid-IR survey do not allow us to properly constrain the 
faint end of the galaxy stellar mass function at such high 
redshifts. This means that a more precise determination of 
the slope a will only be possible once the galaxy stellar mass 
function faint end can be better studied with deeper mid-IR 
surveys. 

• our derived comoving stellar mass densities are (3.0 ± 
1.2) x 10 7 and (2.36ti 8 6 9 8 ) x 10 6 M Mpc" 3 at redshifts z w 
3.25 and 4.6, respectively, which are around 6 and 0.5% of 
the local value. By considering the results of the most recent 
near- and mid-IR surveys, we found that the redshift evo- 
lution of the comoving stellar mass density from z = to 5 
can be modelled as: log 10 p M = -(0.05 ± 0.09) z 2 - (0.22 =F 
0.32) z + 8.69. 



Another key result of our work is the absence of mas- 
sive galaxies at redshifts z > 5. Within our surveyed area of 
0.6 deg 2 , we find only two quite secure candidates at such 
high redshifts, and only one with stellar mass M > 10 11 M©. 
Instead, optical surveys have discovered a substantial num- 
ber of intermediate stellar-mass sources at these redshifts. 
These findings strongly suggest that massive galaxies as a 
significant population only appear at later times, and that 
the epoch around redshifts z ~ 3 — 6 is critical to understand 
the formation of the first massive systems. 

In addition, we have found that a significant fraction of 
the most massive galaxies present at 3 ^ z < 4 would be 
missed by optical surveys, even as deep as R < 27 or z < 26 
magnitudes. This is the consequence of massive galaxies al- 
ready undergoing the main epoch of obscured star formation 
and nuclear activity at these redshifts, which lasts down to 
redshift z ~ 2. 

Globally, this work has shown that deep mid-IR sur- 
veys have a unique role for investigating the first instances 
of massive galaxy assembly at high redshifts, both through 
constraining the rest-frame near-IR galaxy light and discov- 
ering highly reddened objects. Studying significant areas of 
the sky is also key to tracing the most massive galaxies, as 
they quickly become very rare at high redshifts. 
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APPENDIX A: THE IMPACT OF SOURCES 
WITH DEGENERATE SOLUTIONS IN 
REDSHIFT SPACE ON THE GALAXY 
STELLAR MASS FUNCTION AT 3 «: z < 5 

As we explained in Section ^. 21 a total of 244 out of the 1608 
galaxies in our initial z ^ 3 candidate sample have been dis- 
carded on the basis of having degenerate solutions in redshift 
space, i.e. a redshift z < 3 could not be rejected within 3<r 
confidence. The fact that these 244 sources are probably not 
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genuine z > 3 galaxies was supported by the z p hot — z spcc 
diagnostic diagram (see Fig. [2}, as the percentage of outliers 
in it would significantly increase if one considered the high- z 
photometric redshift solutions for these sources. 

In this Appendix our aim is to investigate the impact 
that these discarded sources would have had in our stellar 
mass function at 3 ^ z < 5. 

Fig. lAll shows the stellar mass function that results from 
the STY maximum likelihood analysis when taking into ac- 
count the 244 unlike z 3 candidates along with our final 
sample of genuine z 3 galaxies. 

We see that, for our stellar mass function at 3.0 ^ 
z < 3.5 and 3.5 ^ z < 4.25, the effect of adding these 
extra sources is almost negligible. The values we obtain 
for the Schechter function slope and characteristic stellar 
mass are: a = 1.84±g;°£, M* = (3.16±J;|J) x 10 n M Q , and 
a = 2.10±g;g|, M* = (2.S2toll) x 1O 11 M , in the two 
redshift bins, respectively. These values are similar to those 
quoted in Table [2] This is indicating that the lower redshift 
contaminants have little impact on our stellar mass func- 
tion determinations at 3.0 ^ z < 4.25, and would evenly 
populate the galaxy stellar mass function at different stellar 
masses. 

In the 4.25 ^ z < 5.0, instead, the resulting stellar mass 
function is substantially different to that obtained from the 
genuine z ^ 3 sample. This is because the sample of 244 
unlike z 3 candidates contains 24 sources with 4.25 ^ z < 
5.0, all of which would have corresponding stellar masses 
M > 5 x 10 10 M Q . The addition of these 24 sources triples 
the population of such massive galaxies at 4.25 ^ z < 5.0. 

The best Schechter function obtained through the max- 
imum likelihood analysis on the genuine sample and con- 
taminants altogether yields: a = 2.04j^Q*4 and M* — 
(2.9±?;*) x 10 12 Mq! This M* value is larger than the stellar 
mass of any known galaxy at any redshift. Also, such a stel- 
lar mass function would imply that there are more galaxies 
with stellar masses M > 2 x 10 11 Mq at 4.25 < z < 5.0 than 
at 3.5 ^ z < 4.25. These facts strongly indicate that many 
of the candidates added to the mass function high-mass end 
are very likely not genuine. 
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Figure Al. Comparison of the galaxy stellar mass function ob- 
tained from our secure 2^3 sample (solid lines; see Section 
[4}, with that resulting from incorporating also the extra, unlike 
2^3 candidates (dashed lines). These plots illustrate the effect 
that spurious sources have on the derivation of the stellar mass 
function. 
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